function M = mass_bb(dofs,V,T,d)
%        M = mass_bb(dofs,V,T,d)
% This is normal version, the degree is fixed for all elements
%

% get the areas of all triangles.
x21 = V(T(:,2),1) - V(T(:,1),1);
x31 = V(T(:,3),1) - V(T(:,1),1);
y21 = V(T(:,2),2) - V(T(:,1),2);
y31 = V(T(:,3),2) - V(T(:,1),2);
areas = (x21.*y31 - x31.*y21)/2;

mat = build(d);  [i,j,s] = find(mat);  L = length(i);
m = (d+1)*(d+2)/2;  m = m*m*posT; 
S = zeros(m,1);
Indx1 = zeros(m,1); Indx2 = zeros(m,1); 
pos = 1; % it may be accerate by cancel the loop.
for tri = 1:size(T,1);    
    Indx1(pos:(pos + L-1)) = dofs(tri,i);
    Indx2(pos:(pos + L-1)) = dofs(tri,j);
    S(pos:(pos + L-1)) = areas(tri)*s;
    pos = pos + L;    
end
dim = max(max(dofs));
M = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);